Go to the end of this file for the parameter estimates of the best fitted model on phenotype lobe number
fitting non-linear model with BRMS (Bayesian Regression Models using Stan)
lobe.n <- read_csv("lobe.n.csv", col_types = list(col_character(), col_date(), col_double(), col_double()))
lobe.n
sample.list <- unique(lobe.n$ID)
166 individuals at 4 time points
ggplot(data=lobe.n, aes(x=days, y=lobe_n, group=ID)) +
geom_line(alpha=.1) +
geom_point(size=.1, alpha=.05) +
ggtitle("lobe_number by days")
#### Try Gompertz (4-parameter) Model (see http://www.pisces-conservation.com/growthhelp/index.html?weibul.html) by using BRMS
First set up the formula
gompertz_4p.bf1 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k*(days - I))),
Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
nl=TRUE)
Hmax, asymptotic height at which growth is zero
Hmin, lower asymptotic height
k, growth rate
I, time at the inflection point
stat.data <- function(data) {
data %>%
group_by(days) %>%
summarize(meadian=median(lobe_n),
max=max(lobe_n),
min=min(lobe_n),
sd=sd(lobe_n))
}
stat.data(lobe.n)
Priors. Hmin and Hmax use median at start and end dates.
prior1 <- c(prior(normal(11,4), nlpar="Hmax"),
prior(normal(0,3), nlpar="Hmin"),
prior(normal(1,1), nlpar="k"),
prior(normal(150,10), nlpar="I"))
fit1 <- brm(formula=gompertz_4p.bf1,
data=lobe.n,
prior=prior1)
summary(fit1, waic=TRUE, R2=TRUE)
The model has not converged (some Rhats are > 1.1). Do not analyse the results!
We recommend running more iterations and/or setting stronger priors.There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3598.29; R2 = 0.56
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 10.99 0.46 10.15 11.87 4 1.35
Hmin_Intercept 1.08 0.29 0.53 1.65 122 1.04
k_Intercept 0.29 0.35 0.06 0.93 2 8.46
I_Intercept 136.71 1.53 134.93 140.09 5 1.27
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.63 0.10 3.45 3.82 106 1.03
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit1)
pairs(fit1)
keep k positive. Or more tightly constrain the priors on Hmax and Hmin.
summary(fit2, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3595.49; R2 = 0.56
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 11.17 0.39 10.44 12.01 1218 1.00
Hmin_Intercept 1.07 0.29 0.49 1.62 2396 1.00
k_Intercept 0.09 0.04 0.06 0.21 464 1.01
I_Intercept 137.24 1.37 134.88 140.09 1536 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.61 0.10 3.42 3.82 3255 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit2)
pairs(fit2)
to make k a little bit reasonable
gompertz_4p.bf2 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I))),
Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
nl=TRUE)
summary(fit3, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3595.29; R2 = 0.56
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 11.13 0.37 10.43 11.89 1847 1.00
Hmin_Intercept 1.01 0.30 0.42 1.58 2465 1.00
k_Intercept 0.84 0.23 0.56 1.43 1332 1.00
I_Intercept 136.80 1.30 134.40 139.49 2425 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.61 0.10 3.42 3.81 3699 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit3)
pairs(fit3)
What does the fit look like?
newdata <- data.frame(days=seq(min(lobe.n$days), max(lobe.n$days),1))
fit3.fitted <- cbind(newdata, fitted(fit3, newdata)) %>% as.tibble() %>%
rename(lobe_n=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')
plot
pl <- ggplot(aes(x=days, y=lobe_n),data=NULL)
pl <- pl + geom_line(aes(group=ID), alpha=.1, data=lobe.n)
pl + geom_line(color="skyblue", lwd=1.5, data=fit3.fitted)
Now try adding random effects for model parameters:
What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:
gompertz_4p.fn <- function(Hmax, Hmin, k, I, days) {
Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I)))
}
for(I1 in seq(100,170,10)) {
for(k1 in seq(0,1,.2)) {
tmp.lobe.n <- gompertz_4p.fn(Hmax=11,
Hmin=0,
k=k1,
I=I1,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.lobe.n)
p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
geom_line() + ylim(0,20) + ggtitle(paste0("I=",I1," k=",k1))
print(p)
}
}
First try with only fixing Hmin
gompertz_4p.bf3 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10)*(days - I))),
Hmax + k + I ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I)))
Hmax ~ (1 | ID)
k ~ (1 | ID)
I ~ (1 | ID)
Hmin ~ 1
Data: lobe.n (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 2871.71; R2 = 0.91
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 2.27 0.18 1.94 2.64 1743 1.00
sd(k_Intercept) 0.97 0.20 0.58 1.33 557 1.01
sd(I_Intercept) 17.22 1.14 15.14 19.63 1309 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 11.31 0.23 10.85 11.76 915 1.01
k_Intercept 2.24 0.23 1.83 2.71 4562 1.00
I_Intercept 138.75 1.52 135.79 141.77 1390 1.00
Hmin_Intercept 0.42 0.12 0.19 0.66 10000 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 1.65 0.07 1.52 1.80 1493 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get fitted values
b <- lobe.n %>% ungroup()
fit5.fitted <- cbind(b, fitted(fit5)) %>% as.tibble()
fit5.fitted
plot
plot.fitted(fit5.fitted)
plot fitted vs actual:
plot.fitted.actual <- function(fn) {
r_squared <- summary(lm(Estimate ~ lobe_n, data=fn))$adj.r.squared
r_squared <- round(r_squared, digits = 4)
fn %>%
mutate(days=as.factor(days)) %>%
ggplot(aes(x=lobe_n, y=Estimate, shape=days, color=days)) +
geom_point() +
geom_abline(intercept = 0, slope=1) +
scale_x_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18,20)) +
ggtitle(paste0("R2 = ",r_squared))
}
total_sum_of_squared_residuals <- function(fn) {anova(lm(Estimate ~ lobe_n, data=fn))[2,2]}
plot.fitted.actual(fit5.fitted)
cv
SSR.fit <- total_sum_of_squared_residuals(fit5.fitted)
waic.fit <- round(waic(fit5)$waic, digits=2)
kfoldic.fit <- round(kfold(fit5)$kfoldic, digits=2)
lobe.n.mean <- read_csv("lobe.n.mean.csv")
lobe.n.mean
lapply(1:6, function(i) {
ggplot(aes(x=days, group=ID), data=lobe.n.mean) +
geom_line(aes(y=lobe_n), color="red") + #red: modified data
geom_line(aes(y=lobe_n_raw), color="black") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page =i)
}
)
First set up the formula
gompertz_4p.bf1 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k*(days - I))),
Hmax + Hmin + k + I ~ 1, #general model, paramters do not vary for individuals
nl=TRUE)
stat.data(lobe.n.mean)
stat.data(lobe.n)
Priors. Hmin and Hmax use median at start and end dates.
summary(fit1.mean, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3594.58; R2 = 0.6
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 12.58 0.46 11.77 13.55 1520 1.00
Hmin_Intercept 0.84 0.41 -0.13 1.51 1310 1.00
k_Intercept 0.06 0.01 0.04 0.08 1750 1.00
I_Intercept 136.44 1.62 133.03 139.39 1741 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.61 0.10 3.43 3.82 2901 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit1.mean)
pairs(fit1.mean)
summary(fit2.mean, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-k * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3594.93; R2 = 0.6
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 12.58 0.49 11.75 13.68 1262 1.00
Hmin_Intercept 0.83 0.46 -0.24 1.52 977 1.01
k_Intercept 0.06 0.01 0.04 0.08 1320 1.00
I_Intercept 136.43 1.78 132.72 139.66 1325 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.61 0.10 3.43 3.81 3408 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit2.mean)
pairs(fit2.mean)
gompertz_4p.bf2 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
Hmax + Hmin + k + I ~ 1,
nl=TRUE)
summary(fit3.mean, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I)))
Hmax ~ 1
Hmin ~ 1
k ~ 1
I ~ 1
Data: lobe.n.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = 3594.99; R2 = 0.6
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 12.58 0.48 11.76 13.68 1377 1.00
Hmin_Intercept 0.81 0.47 -0.39 1.52 767 1.00
k_Intercept 0.60 0.11 0.41 0.83 1369 1.00
I_Intercept 136.34 1.82 131.76 139.45 1032 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 3.61 0.10 3.43 3.81 2695 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit3.mean)
pairs(fit3.mean)
What does the fit look like?
newdata <- data.frame(days=seq(min(lobe.n$days), max(lobe.n$days),1))
fit3.mean.fitted <- cbind(newdata, fitted(fit3.mean, newdata)) %>% as.tibble() %>%
rename(lobe_n=Estimate, lower.ci='2.5%ile', upper.ci='97.5%ile')
plot
ggplot(aes(x=days, y=lobe_n),data=NULL) +
geom_line(aes(group=ID), alpha=.1, data=lobe.n) +
geom_line(color="skyblue", lwd=1.5, data=fit3.mean.fitted)
What parameters do we think might be interesting to allow to vary? Probably not Hmin. Try making a series of plots to see how varying delta or k affects things:
lobe.n.fn <- function (Hmax, Hmin, k, delta, days) {
Hmax - (Hmax - Hmin) * exp(-(k/10^6) * (days^delta))
}
for(delta1 in seq(2,4,.25)) {
tmp.lobe.n <- lobe.n.fn(Hmax=8,
Hmin=-2,
k=.29,
delta=delta1,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.lobe.n)
p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
geom_line() + ylim(0,40) + ggtitle("delta=",delta1)
print(p)
}
for(k1 in seq(0,1,.25)) {
tmp.lobe.n <- lobe.n.fn(Hmax=8,
Hmin=-2,
k=k1,
delta=3,
days=newdata$days)
abc <- data.frame(newdata$days, tmp.lobe.n)
p <- ggplot(data=abc, aes(newdata$days, tmp.lobe.n)) +
geom_line() + ylim(0,40) + ggtitle("k=",k1)
print(p)
}
gompertz_4p.bf3 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
Hmax + k + I ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5.mean, waic=TRUE, R2=TRUE)
There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I)))
Hmax ~ (1 | ID)
k ~ (1 | ID)
I ~ (1 | ID)
Hmin ~ 1
Data: lobe.n.mean (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 1803.34; R2 = 0.99
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 2.88 0.17 2.56 3.25 1208 1.01
sd(k_Intercept) 0.94 0.10 0.75 1.16 1959 1.00
sd(I_Intercept) 18.12 1.10 16.06 20.37 1284 1.01
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 12.63 0.24 12.14 13.10 729 1.01
k_Intercept 1.57 0.12 1.35 1.83 2552 1.00
I_Intercept 136.87 1.43 134.11 139.73 465 1.01
Hmin_Intercept 0.06 0.05 -0.04 0.16 10000 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 0.67 0.03 0.61 0.72 2403 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get fitted values
b <- lobe.n.mean %>% ungroup()
fit5.mean.fitted <- cbind(b, fitted(fit5.mean)) %>% as.tibble()
fit5.mean.fitted
plot
plot.fitted(fit5.mean.fitted)
plot fitted vs actual(modified)
plot.fitted.actual(fit5.mean.fitted)
plot fitted vs actual(raw)
plot.fitted.actual.raw(fit5.mean.fitted)
SSR on raw data set
total_sum_of_squared_residuals_raw <- function(fn) {anova(lm(Estimate ~ lobe_n_raw, data=fn))[2,2]}
SSR.fit.mean.raw <- total_sum_of_squared_residuals_raw(fit5.mean.fitted)
cv
SSR.fit.mean <- total_sum_of_squared_residuals(fit5.mean.fitted)
WAIC.fit.mean <- round(waic(fit5.mean)$waic, digits=2)
KFOLDIC.fit.mean <- round(kfold(fit5.mean)$kfoldic, digits=2)
lobe.n.pmm <- read_csv("lobe.n.pmm.csv")
lobe.n.pmm
lapply(1:6, function(i) {
ggplot(aes(x=days, group=ID), data=lobe.n.pmm) +
geom_line(aes(y=lobe_n), color="red") + #red: modified data
geom_line(aes(y=lobe_n_raw), color="black") +
facet_wrap_paginate(~ID, ncol = 6, nrow = 5, page =i)
}
)
stat.data(lobe.n.pmm)
stat.data(lobe.n.mean)
stat.data(lobe.n)
gompertz_4p.bf3 <- bf(
lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I))),
Hmax + k + I ~ (1|ID), # vary for individuals
Hmin ~ 1, # do not vary per individual
nl=TRUE)
summary(fit5.pmm, waic=TRUE, R2=TRUE)
Family: gaussian(identity)
Formula: lobe_n ~ Hmin + (Hmax - Hmin) * exp(-exp(-(k/10) * (days - I)))
Hmax ~ (1 | ID)
k ~ (1 | ID)
I ~ (1 | ID)
Hmin ~ 1
Data: lobe.n.pmm (Number of observations: 664)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
ICs: LOO = NA; WAIC = 2222.45; R2 = 0.97
Group-Level Effects:
~ID (Number of levels: 166)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Hmax_Intercept) 2.78 0.18 2.46 3.15 1954 1.00
sd(k_Intercept) 0.90 0.12 0.69 1.17 2514 1.00
sd(I_Intercept) 17.50 1.07 15.48 19.75 1565 1.01
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Hmax_Intercept 12.50 0.24 12.03 12.97 1614 1.00
k_Intercept 1.80 0.16 1.51 2.14 5011 1.00
I_Intercept 138.42 1.45 135.58 141.29 1015 1.00
Hmin_Intercept 0.22 0.07 0.08 0.36 10000 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 0.97 0.04 0.89 1.06 2727 1.00
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
get fitted values
b <- lobe.n.pmm %>% ungroup()
fit5.pmm.fitted <- cbind(b, fitted(fit5.pmm)) %>% as.tibble()
fit5.pmm.fitted
plot
plot.fitted(fit5.pmm.fitted)
plot fitted vs actual(modified):
plot.fitted.actual(fit5.pmm.fitted)
plot fitted vs actual(raw):
plot.fitted.actual.raw(fit5.pmm.fitted)
SSR on raw data set
SSR.fit.pmm.raw <- total_sum_of_squared_residuals_raw(fit5.pmm.fitted)
cv
SSR.fit.pmm <- total_sum_of_squared_residuals(fit5.pmm.fitted)
WAIC.fit.pmm <- round(waic(fit5.pmm)$waic, digits=2)
KFOLDIC.fit.pmm <- round(kfold(fit5.pmm)$kfoldic, digits=2)
create an empty list for models
cv.list[2,1] <- SSR.fit.mean.raw
cv.list[3,1] <- SSR.fit.pmm.raw
cv.list[1,2] <- SSR.fit
cv.list[2,2] <- SSR.fit.mean
cv.list[3,2] <- SSR.fit.pmm
WAIC(Widely Applicable Information Criterion) an extension of the Akaike Information Criterion (AIC) that is more fully Bayesian.
K-fold CV Data are randomly partitioned into K subsets of equal size. Then the model is refit 10 times(default), each time leaving out one of the 10 subsets.
cv.list[1,3] <- waic.fit
cv.list[2,3] <- WAIC.fit.mean
cv.list[3,3] <- WAIC.fit.pmm
cv.list[1,4] <- kfoldic.fit
cv.list[2,4] <- KFOLDIC.fit.mean
cv.list[3,4] <- KFOLDIC.fit.pmm
cv.list
read
1) population (fixed) & group level (random) effects
2) parameters: Hmax, k & I of each genotype
lobe.n.best.fitted <- fit5.mean$fit
dimnames(lobe.n.best.fitted)
$iterations
NULL
$chains
[1] "chain:1" "chain:2" "chain:3" "chain:4"
$parameters
[1] "b_Hmax_Intercept" "b_k_Intercept" "b_I_Intercept"
[4] "b_Hmin_Intercept" "sd_ID__Hmax_Intercept" "sd_ID__k_Intercept"
[7] "sd_ID__I_Intercept" "sigma" "r_ID__Hmax[ID_1,Intercept]"
[10] "r_ID__Hmax[ID_10,Intercept]" "r_ID__Hmax[ID_100,Intercept]" "r_ID__Hmax[ID_101,Intercept]"
[13] "r_ID__Hmax[ID_102,Intercept]" "r_ID__Hmax[ID_103,Intercept]" "r_ID__Hmax[ID_104,Intercept]"
[16] "r_ID__Hmax[ID_105,Intercept]" "r_ID__Hmax[ID_106,Intercept]" "r_ID__Hmax[ID_107,Intercept]"
[19] "r_ID__Hmax[ID_109,Intercept]" "r_ID__Hmax[ID_111,Intercept]" "r_ID__Hmax[ID_112,Intercept]"
[22] "r_ID__Hmax[ID_113,Intercept]" "r_ID__Hmax[ID_114,Intercept]" "r_ID__Hmax[ID_116,Intercept]"
[25] "r_ID__Hmax[ID_117,Intercept]" "r_ID__Hmax[ID_118,Intercept]" "r_ID__Hmax[ID_119,Intercept]"
[28] "r_ID__Hmax[ID_12,Intercept]" "r_ID__Hmax[ID_120,Intercept]" "r_ID__Hmax[ID_121,Intercept]"
[31] "r_ID__Hmax[ID_122,Intercept]" "r_ID__Hmax[ID_123,Intercept]" "r_ID__Hmax[ID_124,Intercept]"
[34] "r_ID__Hmax[ID_125,Intercept]" "r_ID__Hmax[ID_126,Intercept]" "r_ID__Hmax[ID_127,Intercept]"
[37] "r_ID__Hmax[ID_128,Intercept]" "r_ID__Hmax[ID_129,Intercept]" "r_ID__Hmax[ID_13,Intercept]"
[40] "r_ID__Hmax[ID_130,Intercept]" "r_ID__Hmax[ID_132,Intercept]" "r_ID__Hmax[ID_133,Intercept]"
[43] "r_ID__Hmax[ID_134,Intercept]" "r_ID__Hmax[ID_135,Intercept]" "r_ID__Hmax[ID_136,Intercept]"
[46] "r_ID__Hmax[ID_137,Intercept]" "r_ID__Hmax[ID_138,Intercept]" "r_ID__Hmax[ID_139,Intercept]"
[49] "r_ID__Hmax[ID_14,Intercept]" "r_ID__Hmax[ID_142,Intercept]" "r_ID__Hmax[ID_143,Intercept]"
[52] "r_ID__Hmax[ID_144,Intercept]" "r_ID__Hmax[ID_147,Intercept]" "r_ID__Hmax[ID_148,Intercept]"
[55] "r_ID__Hmax[ID_15,Intercept]" "r_ID__Hmax[ID_152,Intercept]" "r_ID__Hmax[ID_153,Intercept]"
[58] "r_ID__Hmax[ID_154,Intercept]" "r_ID__Hmax[ID_155,Intercept]" "r_ID__Hmax[ID_156,Intercept]"
[61] "r_ID__Hmax[ID_157,Intercept]" "r_ID__Hmax[ID_158,Intercept]" "r_ID__Hmax[ID_16,Intercept]"
[64] "r_ID__Hmax[ID_160,Intercept]" "r_ID__Hmax[ID_161,Intercept]" "r_ID__Hmax[ID_163,Intercept]"
[67] "r_ID__Hmax[ID_165,Intercept]" "r_ID__Hmax[ID_166,Intercept]" "r_ID__Hmax[ID_169,Intercept]"
[70] "r_ID__Hmax[ID_17,Intercept]" "r_ID__Hmax[ID_170,Intercept]" "r_ID__Hmax[ID_171,Intercept]"
[73] "r_ID__Hmax[ID_172,Intercept]" "r_ID__Hmax[ID_174,Intercept]" "r_ID__Hmax[ID_175,Intercept]"
[76] "r_ID__Hmax[ID_176,Intercept]" "r_ID__Hmax[ID_177,Intercept]" "r_ID__Hmax[ID_18,Intercept]"
[79] "r_ID__Hmax[ID_180,Intercept]" "r_ID__Hmax[ID_181,Intercept]" "r_ID__Hmax[ID_186,Intercept]"
[82] "r_ID__Hmax[ID_189,Intercept]" "r_ID__Hmax[ID_19,Intercept]" "r_ID__Hmax[ID_191,Intercept]"
[85] "r_ID__Hmax[ID_192,Intercept]" "r_ID__Hmax[ID_197,Intercept]" "r_ID__Hmax[ID_198,Intercept]"
[88] "r_ID__Hmax[ID_199,Intercept]" "r_ID__Hmax[ID_2,Intercept]" "r_ID__Hmax[ID_20,Intercept]"
[91] "r_ID__Hmax[ID_201,Intercept]" "r_ID__Hmax[ID_202,Intercept]" "r_ID__Hmax[ID_207,Intercept]"
[94] "r_ID__Hmax[ID_208,Intercept]" "r_ID__Hmax[ID_21,Intercept]" "r_ID__Hmax[ID_212,Intercept]"
[97] "r_ID__Hmax[ID_214,Intercept]" "r_ID__Hmax[ID_217,Intercept]" "r_ID__Hmax[ID_223,Intercept]"
[100] "r_ID__Hmax[ID_224,Intercept]" "r_ID__Hmax[ID_227,Intercept]" "r_ID__Hmax[ID_23,Intercept]"
[103] "r_ID__Hmax[ID_231,Intercept]" "r_ID__Hmax[ID_234,Intercept]" "r_ID__Hmax[ID_24,Intercept]"
[106] "r_ID__Hmax[ID_25,Intercept]" "r_ID__Hmax[ID_26,Intercept]" "r_ID__Hmax[ID_27,Intercept]"
[109] "r_ID__Hmax[ID_28,Intercept]" "r_ID__Hmax[ID_29,Intercept]" "r_ID__Hmax[ID_3,Intercept]"
[112] "r_ID__Hmax[ID_30,Intercept]" "r_ID__Hmax[ID_32,Intercept]" "r_ID__Hmax[ID_34,Intercept]"
[115] "r_ID__Hmax[ID_35,Intercept]" "r_ID__Hmax[ID_36,Intercept]" "r_ID__Hmax[ID_37,Intercept]"
[118] "r_ID__Hmax[ID_38,Intercept]" "r_ID__Hmax[ID_39,Intercept]" "r_ID__Hmax[ID_4,Intercept]"
[121] "r_ID__Hmax[ID_40,Intercept]" "r_ID__Hmax[ID_41,Intercept]" "r_ID__Hmax[ID_42,Intercept]"
[124] "r_ID__Hmax[ID_45,Intercept]" "r_ID__Hmax[ID_47,Intercept]" "r_ID__Hmax[ID_48,Intercept]"
[127] "r_ID__Hmax[ID_49,Intercept]" "r_ID__Hmax[ID_50,Intercept]" "r_ID__Hmax[ID_51,Intercept]"
[130] "r_ID__Hmax[ID_52,Intercept]" "r_ID__Hmax[ID_53,Intercept]" "r_ID__Hmax[ID_55,Intercept]"
[133] "r_ID__Hmax[ID_56,Intercept]" "r_ID__Hmax[ID_57,Intercept]" "r_ID__Hmax[ID_58,Intercept]"
[136] "r_ID__Hmax[ID_59,Intercept]" "r_ID__Hmax[ID_6,Intercept]" "r_ID__Hmax[ID_60,Intercept]"
[139] "r_ID__Hmax[ID_61,Intercept]" "r_ID__Hmax[ID_62,Intercept]" "r_ID__Hmax[ID_63,Intercept]"
[142] "r_ID__Hmax[ID_64,Intercept]" "r_ID__Hmax[ID_65,Intercept]" "r_ID__Hmax[ID_67,Intercept]"
[145] "r_ID__Hmax[ID_68,Intercept]" "r_ID__Hmax[ID_69,Intercept]" "r_ID__Hmax[ID_7,Intercept]"
[148] "r_ID__Hmax[ID_72,Intercept]" "r_ID__Hmax[ID_73,Intercept]" "r_ID__Hmax[ID_74,Intercept]"
[151] "r_ID__Hmax[ID_75,Intercept]" "r_ID__Hmax[ID_76,Intercept]" "r_ID__Hmax[ID_77,Intercept]"
[154] "r_ID__Hmax[ID_78,Intercept]" "r_ID__Hmax[ID_79,Intercept]" "r_ID__Hmax[ID_8,Intercept]"
[157] "r_ID__Hmax[ID_80,Intercept]" "r_ID__Hmax[ID_81,Intercept]" "r_ID__Hmax[ID_82,Intercept]"
[160] "r_ID__Hmax[ID_84,Intercept]" "r_ID__Hmax[ID_85,Intercept]" "r_ID__Hmax[ID_86,Intercept]"
[163] "r_ID__Hmax[ID_87,Intercept]" "r_ID__Hmax[ID_88,Intercept]" "r_ID__Hmax[ID_89,Intercept]"
[166] "r_ID__Hmax[ID_9,Intercept]" "r_ID__Hmax[ID_90,Intercept]" "r_ID__Hmax[ID_91,Intercept]"
[169] "r_ID__Hmax[ID_92,Intercept]" "r_ID__Hmax[ID_94,Intercept]" "r_ID__Hmax[ID_95,Intercept]"
[172] "r_ID__Hmax[ID_96,Intercept]" "r_ID__Hmax[ID_98,Intercept]" "r_ID__Hmax[ID_99,Intercept]"
[175] "r_ID__k[ID_1,Intercept]" "r_ID__k[ID_10,Intercept]" "r_ID__k[ID_100,Intercept]"
[178] "r_ID__k[ID_101,Intercept]" "r_ID__k[ID_102,Intercept]" "r_ID__k[ID_103,Intercept]"
[181] "r_ID__k[ID_104,Intercept]" "r_ID__k[ID_105,Intercept]" "r_ID__k[ID_106,Intercept]"
[184] "r_ID__k[ID_107,Intercept]" "r_ID__k[ID_109,Intercept]" "r_ID__k[ID_111,Intercept]"
[187] "r_ID__k[ID_112,Intercept]" "r_ID__k[ID_113,Intercept]" "r_ID__k[ID_114,Intercept]"
[190] "r_ID__k[ID_116,Intercept]" "r_ID__k[ID_117,Intercept]" "r_ID__k[ID_118,Intercept]"
[193] "r_ID__k[ID_119,Intercept]" "r_ID__k[ID_12,Intercept]" "r_ID__k[ID_120,Intercept]"
[196] "r_ID__k[ID_121,Intercept]" "r_ID__k[ID_122,Intercept]" "r_ID__k[ID_123,Intercept]"
[199] "r_ID__k[ID_124,Intercept]" "r_ID__k[ID_125,Intercept]" "r_ID__k[ID_126,Intercept]"
[202] "r_ID__k[ID_127,Intercept]" "r_ID__k[ID_128,Intercept]" "r_ID__k[ID_129,Intercept]"
[205] "r_ID__k[ID_13,Intercept]" "r_ID__k[ID_130,Intercept]" "r_ID__k[ID_132,Intercept]"
[208] "r_ID__k[ID_133,Intercept]" "r_ID__k[ID_134,Intercept]" "r_ID__k[ID_135,Intercept]"
[211] "r_ID__k[ID_136,Intercept]" "r_ID__k[ID_137,Intercept]" "r_ID__k[ID_138,Intercept]"
[214] "r_ID__k[ID_139,Intercept]" "r_ID__k[ID_14,Intercept]" "r_ID__k[ID_142,Intercept]"
[217] "r_ID__k[ID_143,Intercept]" "r_ID__k[ID_144,Intercept]" "r_ID__k[ID_147,Intercept]"
[220] "r_ID__k[ID_148,Intercept]" "r_ID__k[ID_15,Intercept]" "r_ID__k[ID_152,Intercept]"
[223] "r_ID__k[ID_153,Intercept]" "r_ID__k[ID_154,Intercept]" "r_ID__k[ID_155,Intercept]"
[226] "r_ID__k[ID_156,Intercept]" "r_ID__k[ID_157,Intercept]" "r_ID__k[ID_158,Intercept]"
[229] "r_ID__k[ID_16,Intercept]" "r_ID__k[ID_160,Intercept]" "r_ID__k[ID_161,Intercept]"
[232] "r_ID__k[ID_163,Intercept]" "r_ID__k[ID_165,Intercept]" "r_ID__k[ID_166,Intercept]"
[235] "r_ID__k[ID_169,Intercept]" "r_ID__k[ID_17,Intercept]" "r_ID__k[ID_170,Intercept]"
[238] "r_ID__k[ID_171,Intercept]" "r_ID__k[ID_172,Intercept]" "r_ID__k[ID_174,Intercept]"
[241] "r_ID__k[ID_175,Intercept]" "r_ID__k[ID_176,Intercept]" "r_ID__k[ID_177,Intercept]"
[244] "r_ID__k[ID_18,Intercept]" "r_ID__k[ID_180,Intercept]" "r_ID__k[ID_181,Intercept]"
[247] "r_ID__k[ID_186,Intercept]" "r_ID__k[ID_189,Intercept]" "r_ID__k[ID_19,Intercept]"
[250] "r_ID__k[ID_191,Intercept]" "r_ID__k[ID_192,Intercept]" "r_ID__k[ID_197,Intercept]"
[253] "r_ID__k[ID_198,Intercept]" "r_ID__k[ID_199,Intercept]" "r_ID__k[ID_2,Intercept]"
[256] "r_ID__k[ID_20,Intercept]" "r_ID__k[ID_201,Intercept]" "r_ID__k[ID_202,Intercept]"
[259] "r_ID__k[ID_207,Intercept]" "r_ID__k[ID_208,Intercept]" "r_ID__k[ID_21,Intercept]"
[262] "r_ID__k[ID_212,Intercept]" "r_ID__k[ID_214,Intercept]" "r_ID__k[ID_217,Intercept]"
[265] "r_ID__k[ID_223,Intercept]" "r_ID__k[ID_224,Intercept]" "r_ID__k[ID_227,Intercept]"
[268] "r_ID__k[ID_23,Intercept]" "r_ID__k[ID_231,Intercept]" "r_ID__k[ID_234,Intercept]"
[271] "r_ID__k[ID_24,Intercept]" "r_ID__k[ID_25,Intercept]" "r_ID__k[ID_26,Intercept]"
[274] "r_ID__k[ID_27,Intercept]" "r_ID__k[ID_28,Intercept]" "r_ID__k[ID_29,Intercept]"
[277] "r_ID__k[ID_3,Intercept]" "r_ID__k[ID_30,Intercept]" "r_ID__k[ID_32,Intercept]"
[280] "r_ID__k[ID_34,Intercept]" "r_ID__k[ID_35,Intercept]" "r_ID__k[ID_36,Intercept]"
[283] "r_ID__k[ID_37,Intercept]" "r_ID__k[ID_38,Intercept]" "r_ID__k[ID_39,Intercept]"
[286] "r_ID__k[ID_4,Intercept]" "r_ID__k[ID_40,Intercept]" "r_ID__k[ID_41,Intercept]"
[289] "r_ID__k[ID_42,Intercept]" "r_ID__k[ID_45,Intercept]" "r_ID__k[ID_47,Intercept]"
[292] "r_ID__k[ID_48,Intercept]" "r_ID__k[ID_49,Intercept]" "r_ID__k[ID_50,Intercept]"
[295] "r_ID__k[ID_51,Intercept]" "r_ID__k[ID_52,Intercept]" "r_ID__k[ID_53,Intercept]"
[298] "r_ID__k[ID_55,Intercept]" "r_ID__k[ID_56,Intercept]" "r_ID__k[ID_57,Intercept]"
[301] "r_ID__k[ID_58,Intercept]" "r_ID__k[ID_59,Intercept]" "r_ID__k[ID_6,Intercept]"
[304] "r_ID__k[ID_60,Intercept]" "r_ID__k[ID_61,Intercept]" "r_ID__k[ID_62,Intercept]"
[307] "r_ID__k[ID_63,Intercept]" "r_ID__k[ID_64,Intercept]" "r_ID__k[ID_65,Intercept]"
[310] "r_ID__k[ID_67,Intercept]" "r_ID__k[ID_68,Intercept]" "r_ID__k[ID_69,Intercept]"
[313] "r_ID__k[ID_7,Intercept]" "r_ID__k[ID_72,Intercept]" "r_ID__k[ID_73,Intercept]"
[316] "r_ID__k[ID_74,Intercept]" "r_ID__k[ID_75,Intercept]" "r_ID__k[ID_76,Intercept]"
[319] "r_ID__k[ID_77,Intercept]" "r_ID__k[ID_78,Intercept]" "r_ID__k[ID_79,Intercept]"
[322] "r_ID__k[ID_8,Intercept]" "r_ID__k[ID_80,Intercept]" "r_ID__k[ID_81,Intercept]"
[325] "r_ID__k[ID_82,Intercept]" "r_ID__k[ID_84,Intercept]" "r_ID__k[ID_85,Intercept]"
[328] "r_ID__k[ID_86,Intercept]" "r_ID__k[ID_87,Intercept]" "r_ID__k[ID_88,Intercept]"
[331] "r_ID__k[ID_89,Intercept]" "r_ID__k[ID_9,Intercept]" "r_ID__k[ID_90,Intercept]"
[334] "r_ID__k[ID_91,Intercept]" "r_ID__k[ID_92,Intercept]" "r_ID__k[ID_94,Intercept]"
[337] "r_ID__k[ID_95,Intercept]" "r_ID__k[ID_96,Intercept]" "r_ID__k[ID_98,Intercept]"
[340] "r_ID__k[ID_99,Intercept]" "r_ID__I[ID_1,Intercept]" "r_ID__I[ID_10,Intercept]"
[343] "r_ID__I[ID_100,Intercept]" "r_ID__I[ID_101,Intercept]" "r_ID__I[ID_102,Intercept]"
[346] "r_ID__I[ID_103,Intercept]" "r_ID__I[ID_104,Intercept]" "r_ID__I[ID_105,Intercept]"
[349] "r_ID__I[ID_106,Intercept]" "r_ID__I[ID_107,Intercept]" "r_ID__I[ID_109,Intercept]"
[352] "r_ID__I[ID_111,Intercept]" "r_ID__I[ID_112,Intercept]" "r_ID__I[ID_113,Intercept]"
[355] "r_ID__I[ID_114,Intercept]" "r_ID__I[ID_116,Intercept]" "r_ID__I[ID_117,Intercept]"
[358] "r_ID__I[ID_118,Intercept]" "r_ID__I[ID_119,Intercept]" "r_ID__I[ID_12,Intercept]"
[361] "r_ID__I[ID_120,Intercept]" "r_ID__I[ID_121,Intercept]" "r_ID__I[ID_122,Intercept]"
[364] "r_ID__I[ID_123,Intercept]" "r_ID__I[ID_124,Intercept]" "r_ID__I[ID_125,Intercept]"
[367] "r_ID__I[ID_126,Intercept]" "r_ID__I[ID_127,Intercept]" "r_ID__I[ID_128,Intercept]"
[370] "r_ID__I[ID_129,Intercept]" "r_ID__I[ID_13,Intercept]" "r_ID__I[ID_130,Intercept]"
[373] "r_ID__I[ID_132,Intercept]" "r_ID__I[ID_133,Intercept]" "r_ID__I[ID_134,Intercept]"
[376] "r_ID__I[ID_135,Intercept]" "r_ID__I[ID_136,Intercept]" "r_ID__I[ID_137,Intercept]"
[379] "r_ID__I[ID_138,Intercept]" "r_ID__I[ID_139,Intercept]" "r_ID__I[ID_14,Intercept]"
[382] "r_ID__I[ID_142,Intercept]" "r_ID__I[ID_143,Intercept]" "r_ID__I[ID_144,Intercept]"
[385] "r_ID__I[ID_147,Intercept]" "r_ID__I[ID_148,Intercept]" "r_ID__I[ID_15,Intercept]"
[388] "r_ID__I[ID_152,Intercept]" "r_ID__I[ID_153,Intercept]" "r_ID__I[ID_154,Intercept]"
[391] "r_ID__I[ID_155,Intercept]" "r_ID__I[ID_156,Intercept]" "r_ID__I[ID_157,Intercept]"
[394] "r_ID__I[ID_158,Intercept]" "r_ID__I[ID_16,Intercept]" "r_ID__I[ID_160,Intercept]"
[397] "r_ID__I[ID_161,Intercept]" "r_ID__I[ID_163,Intercept]" "r_ID__I[ID_165,Intercept]"
[400] "r_ID__I[ID_166,Intercept]" "r_ID__I[ID_169,Intercept]" "r_ID__I[ID_17,Intercept]"
[403] "r_ID__I[ID_170,Intercept]" "r_ID__I[ID_171,Intercept]" "r_ID__I[ID_172,Intercept]"
[406] "r_ID__I[ID_174,Intercept]" "r_ID__I[ID_175,Intercept]" "r_ID__I[ID_176,Intercept]"
[409] "r_ID__I[ID_177,Intercept]" "r_ID__I[ID_18,Intercept]" "r_ID__I[ID_180,Intercept]"
[412] "r_ID__I[ID_181,Intercept]" "r_ID__I[ID_186,Intercept]" "r_ID__I[ID_189,Intercept]"
[415] "r_ID__I[ID_19,Intercept]" "r_ID__I[ID_191,Intercept]" "r_ID__I[ID_192,Intercept]"
[418] "r_ID__I[ID_197,Intercept]" "r_ID__I[ID_198,Intercept]" "r_ID__I[ID_199,Intercept]"
[421] "r_ID__I[ID_2,Intercept]" "r_ID__I[ID_20,Intercept]" "r_ID__I[ID_201,Intercept]"
[424] "r_ID__I[ID_202,Intercept]" "r_ID__I[ID_207,Intercept]" "r_ID__I[ID_208,Intercept]"
[427] "r_ID__I[ID_21,Intercept]" "r_ID__I[ID_212,Intercept]" "r_ID__I[ID_214,Intercept]"
[430] "r_ID__I[ID_217,Intercept]" "r_ID__I[ID_223,Intercept]" "r_ID__I[ID_224,Intercept]"
[433] "r_ID__I[ID_227,Intercept]" "r_ID__I[ID_23,Intercept]" "r_ID__I[ID_231,Intercept]"
[436] "r_ID__I[ID_234,Intercept]" "r_ID__I[ID_24,Intercept]" "r_ID__I[ID_25,Intercept]"
[439] "r_ID__I[ID_26,Intercept]" "r_ID__I[ID_27,Intercept]" "r_ID__I[ID_28,Intercept]"
[442] "r_ID__I[ID_29,Intercept]" "r_ID__I[ID_3,Intercept]" "r_ID__I[ID_30,Intercept]"
[445] "r_ID__I[ID_32,Intercept]" "r_ID__I[ID_34,Intercept]" "r_ID__I[ID_35,Intercept]"
[448] "r_ID__I[ID_36,Intercept]" "r_ID__I[ID_37,Intercept]" "r_ID__I[ID_38,Intercept]"
[451] "r_ID__I[ID_39,Intercept]" "r_ID__I[ID_4,Intercept]" "r_ID__I[ID_40,Intercept]"
[454] "r_ID__I[ID_41,Intercept]" "r_ID__I[ID_42,Intercept]" "r_ID__I[ID_45,Intercept]"
[457] "r_ID__I[ID_47,Intercept]" "r_ID__I[ID_48,Intercept]" "r_ID__I[ID_49,Intercept]"
[460] "r_ID__I[ID_50,Intercept]" "r_ID__I[ID_51,Intercept]" "r_ID__I[ID_52,Intercept]"
[463] "r_ID__I[ID_53,Intercept]" "r_ID__I[ID_55,Intercept]" "r_ID__I[ID_56,Intercept]"
[466] "r_ID__I[ID_57,Intercept]" "r_ID__I[ID_58,Intercept]" "r_ID__I[ID_59,Intercept]"
[469] "r_ID__I[ID_6,Intercept]" "r_ID__I[ID_60,Intercept]" "r_ID__I[ID_61,Intercept]"
[472] "r_ID__I[ID_62,Intercept]" "r_ID__I[ID_63,Intercept]" "r_ID__I[ID_64,Intercept]"
[475] "r_ID__I[ID_65,Intercept]" "r_ID__I[ID_67,Intercept]" "r_ID__I[ID_68,Intercept]"
[478] "r_ID__I[ID_69,Intercept]" "r_ID__I[ID_7,Intercept]" "r_ID__I[ID_72,Intercept]"
[481] "r_ID__I[ID_73,Intercept]" "r_ID__I[ID_74,Intercept]" "r_ID__I[ID_75,Intercept]"
[484] "r_ID__I[ID_76,Intercept]" "r_ID__I[ID_77,Intercept]" "r_ID__I[ID_78,Intercept]"
[487] "r_ID__I[ID_79,Intercept]" "r_ID__I[ID_8,Intercept]" "r_ID__I[ID_80,Intercept]"
[490] "r_ID__I[ID_81,Intercept]" "r_ID__I[ID_82,Intercept]" "r_ID__I[ID_84,Intercept]"
[493] "r_ID__I[ID_85,Intercept]" "r_ID__I[ID_86,Intercept]" "r_ID__I[ID_87,Intercept]"
[496] "r_ID__I[ID_88,Intercept]" "r_ID__I[ID_89,Intercept]" "r_ID__I[ID_9,Intercept]"
[499] "r_ID__I[ID_90,Intercept]" "r_ID__I[ID_91,Intercept]" "r_ID__I[ID_92,Intercept]"
[502] "r_ID__I[ID_94,Intercept]" "r_ID__I[ID_95,Intercept]" "r_ID__I[ID_96,Intercept]"
[505] "r_ID__I[ID_98,Intercept]" "r_ID__I[ID_99,Intercept]" "lp__"
lobe.n.best.fitted.summary <- summary(lobe.n.best.fitted)$summary
write.csv(lobe.n.best.fitted.summary, file = "/Users/seungmokim/414_Growth_Model/summary/growth_model_phenotypes/lobe_number/lobe.n.best.fitted.summary.csv")
r_ID__Hmax
datalist0 = list()
for (samples in sample.list) {
ttt0 <- paste0("r_ID__Hmax[",samples,",Intercept]")
dat0 <- extract(lobe.n.best.fitted, ttt0, permuted=TRUE)
datalist0[[ttt0]] <- dat0
}
list.data0 = do.call(rbind, datalist0)
Hmax.ID <- data.frame(matrix(unlist(list.data0), nrow = 166, byrow = T))
row.names(Hmax.ID) <- paste0("r_ID__Hmax[",sample.list,",Intercept]")
head(Hmax.ID)
r_ID__k
datalist1 = list()
for (samples in sample.list) {
ttt1 <- paste0("r_ID__k[",samples,",Intercept]")
dat1 <- extract(lobe.n.best.fitted, ttt1, permuted=TRUE)
datalist1[[ttt1]] <- dat1
}
list.data1 = do.call(rbind, datalist1)
k.ID <- data.frame(matrix(unlist(list.data1), nrow = 166, byrow = T))
row.names(k.ID) <- paste0("r_ID__k[",sample.list,",Intercept]")
head(k.ID)
r_ID__I
datalist2 = list()
for (samples in sample.list) {
ttt2 <- paste0("r_ID__I[",samples,",Intercept]")
dat2 <- extract(lobe.n.best.fitted, ttt2, permuted=TRUE)
datalist2[[ttt2]] <- dat2
}
list.data2 = do.call(rbind, datalist2)
I.ID <- data.frame(matrix(unlist(list.data2), nrow = 166, byrow = T))
row.names(I.ID) <- paste0("r_ID__I[",sample.list,",Intercept]")
head(I.ID)
estimates of parameters
par.list = c("b_Hmax_Intercept", "b_k_Intercept", "b_I_Intercept", "b_Hmin_Intercept", "sd_ID__Hmax_Intercept", "sd_ID__k_Intercept","sd_ID__I_Intercept","sigma", "lp__")
datalist3 = list()
for (pars in par.list) {
dat3 <- extract(lobe.n.best.fitted, pars, permuted=TRUE)
datalist3[[pars]] <- dat3
}
list.data3 = do.call(rbind, datalist3)
pars.est <- data.frame(matrix(unlist(list.data3), nrow = 9, byrow = T))
row.names(pars.est) <- par.list
pars.est
combine estimates of all parameters
pars.estimate.lobe.n.fitted <- rbind(pars.est, Hmax.ID, k.ID, I.ID)
head(pars.estimate.lobe.n.fitted, 10)
write.csv(pars.estimate.lobe.n.fitted, file="pars.estimate.lobe.n.fitted.csv")